*Program name:      eventtimes_v*.do
*Written by:        Diane Alexander
*Date written:      April 3, 2020
*Last modified by:  Ezra Karger
*Date modified:     July 9, 2021

clear all

*Set sub-directory paths--same for all users:

global unacast  "${machine}/data/clean/unacast"
global stay_at_home "${machine}/data/clean/stay_at_home"
global womply "${machine}/data/clean/womply"
global nyt_cases "${machine}/data/clean/nyt_cases"
global census "${machine}/data/clean/census"
global trump "${machine}/data/clean/trump"
global acs "${machine}/data/clean/acs_2018"
global usda "${machine}/data/clean/usda_ers"
global weather "${machine}/data/clean/weather"
global covid_tracker "${machine}/data/clean/covid_tracker"
global safegraph "${machine}/data/clean/safegraph"

global outdata "${machine}/data/clean/"

global logs "${machine}/logs"
global tables "${machine}/tables/event_time/main"
global figures "${machine}/figures/event_time/main"

capture log close
log using "$eventtimes_no_stateday_fes.log", text replace

*Read in coded stay-at-home orders

use "${stay_at_home}/stay_at_home_AM.dta"


**Merge on Unacast data

merge 1:1 county_fips date using "${unacast}/social_distance_county_04_29.dta"
drop _merge

rename daily_distance_diff distance
rename daily_visitation_diff visit 

**drop dates after April 17 (10 days after last order on April 7)
drop if date > mdy(4,17,2020)


**Merge on Womply data, but only after making Womply's data wide by category
*	of business. The categories include arts and entertainment (like movie theaters),
*	 bars, restaurants, and auto service shops.
*Final Womply data is at the date*county level, with variables for revenue from
*	each business type.
*For this analysis, we combine all businesses into two categories: food-related
*	establishments and other establishments

preserve
use "${womply}/womply_daily_revenue", clear
keep if year(date)==2020
keep if !missing(county_fips)

tab category, m

gen clean_cat="other"
replace clean_cat="food" if inlist(category,"bars and lounges","food and beverage shops","quick serve food and beverage businesses","restaurants")
tab category clean_cat, m


keep date county_fips revenue clean_cat intra_covid_active_merchants intra_covid_inactive_merchants merchants
rename revenue revenue_
rename intra_covid_active_merchants active_
rename intra_covid_inactive_merchants inactive_
rename merchants merchants_
collapse (sum) revenue_ active_ inactive_ merchants_, by(date county clean_cat)

reshape wide revenue_ active_ inactive_ merchants_, i(date county_fips) j(clean_cat) string

*Make 'all revenue' category
egen revenue_all = rowtotal(revenue_*)
egen inactive_all = rowtotal(inactive_*)
egen active_all = rowtotal(active_*)
egen merchants_all = rowtotal(merchants_*)

gen fr_i_all=inactive_all/(merchants_all)
gen fr_i_food=inactive_food/(merchants_food)
gen fr_i_other=inactive_other/(merchants_other)


*Missing observations in a date*county*category are instances with zero
*	revenue for that category of business.

foreach var of varlist revenue* {
	replace `var'=0 if missing(`var')
	
	gen log`var'=log(`var'+1)
}

**drop dates after April 17 (10 days after last order on April 7)
drop if date> mdy(4,17,2020)

tempfile womply
save "`womply'"
restore

merge 1:1 date county_fips using "`womply'", nogen



*Merge on USDA county-level characteristics
**CHECK THIS MERGE
merge m:1 county_fips using "${usda}/county_information.dta"
reg tot_pop pop_estimate_2018
*drop puerto rico
drop if county_fips>57000


*Create revenue per-capita estimates from the Womply panel

gen firstweekmarch=1 if date >= mdy(3,1,2020) & date <= mdy(3,7,2020)

tab date if firstweekmarch==1


foreach v in food all other {
display "`var'"

	gen temp_estab_`v'=(merchants_`v') if firstweekmarch==1
	bysort county_fips: egen estab_earlymarch_`v'=mean(temp_estab_`v')
	gen perest_revenue_`v'=revenue_`v'/(estab_earlymarch_`v')
	drop temp_estab_`v'
}

rename perest_revenue_food pe_r_food
rename perest_revenue_other pe_r_other
rename perest_revenue_all pe_r_all

gen log_pe_food = log(pe_r_food+1)
gen log_pe_other = log(pe_r_other+1)
gen log_pe_all = log(pe_r_all+1)


sort county_fips date
 
*Merge on Trump vote share

preserve
import delimited using "${trump}/2016_US_County_Level_Presidential_Results.csv", clear

rename combined_fips county_fips
format county_fips %05.0f
gen gop_win_pct = per_gop-per_dem
keep county_fips gop_win_pct

sum gop_win_pct, d

tempfile trump
save "`trump'"
restore

merge m:1 county_fips using "`trump'", nogen

gen split_trump=gop_win_pct>0 & !missing(gop_win_pct)
gen split_dem=gop_win_pct<=0 & !missing(gop_win_pct)




**This just pulls the effective date of the stay-at-home order through the whole time period
foreach x in effective_date {
	bysort county_fips: egen `x'_t=mean(`x')
	drop `x'
	rename `x'_t `x'
}
format effective_date %td


** 88 counties don't show up in the social distance data, but these counties have very small populations

bysort county_fips: egen ever_stayathome=max(stay_at_home)

gen weekend=1 if weekday==1 | weekday==7

drop state_fips
gen temp=county_fips
format temp %05.0f
tostring temp,  replace usedisplayformat
gen state_fips=substr(temp ,1,2)
destring state_fips, replace
drop temp state




*Merge on weather data
preserve
use "${weather}/state_daily_2020.dta", clear
desc

rename date_mdy date

sum PRCP, d
gen log_precipitation = log(PRCP+1)
sum log_precipitation, d

keep state_fips date TEMP log_precipitation
keep if inrange(state_fips,1,56)

tempfile weather
save "`weather'"
restore

capture drop _merge
merge m:1 state_fips date using "`weather'"
drop if _merge==2



*Merge on state-date Covid-19 data
preserve
use "${covid_tracker}/covid_tracker_states_5_28.dta", clear
desc

rename date_mdy date

keep state_fips date positive negative hospitalized
gen log_positive = log(positive+1 )
gen log_negative = log(negative+1 )
gen log_hospitalized = log(hospitalized+1 )

keep if inrange(state_fips,1,56) & inrange(date,mdy(3,1,2020),mdy(5,30,2020))

duplicates report state_fips date

tab state_fips, m
tab date, m

tempfile tracker
save "`tracker'"
restore

capture drop _merge
merge m:1 state_fips date using "`tracker'"

tab date _merge, m
drop if _merge==2



**unacast data starts feb 24==21969, but NYTimes data starts in January
merge 1:1 county_fips date using "${nyt_cases}/nyt_cases_04_17.dta", nogen 

**missings in the cases/deaths are zeroes
replace cases=0 if cases==.
replace cases_daily_incr=0 if cases_daily_incr==. | cases_daily_incr<0
replace deaths=0 if deaths==.
replace deaths_daily_incr=0 if deaths_daily_incr==. | deaths_daily_incr<0

*Try controlling for cases and deaths (logged)?
gen log_deaths =  log(deaths_daily_incr+1)
gen log_cases =  log(cases_daily_incr+1)


*Create control variable which is 'days since first case', topcoded at 30.
*Add 100 to the control so that we can include it in regressions as a factor
*	variable, with all positive values.

sort county_fips date
gen case_occurred = date if cases>0
by county_fips: egen date_of_first_case=min(case_occurred)
gen days_since_first_case = date-date_of_first_case
replace days_since_first_case=999 if missing(days_since_first_case)
replace days_since_first_case=30 if inrange(days_since_first_case,30,998)
replace days_since_first_case = days_since_first_case + 100


*Create low- and high-income bins using top and bottom quartile
*	median household income in 2018

preserve
keep county_fips median_household_income_2018
duplicates drop
sum median_household_income_2018, d
drop if median_household_income_2018==.


xtile quantile_inc=median_household_income_2018, n(2)
tab quantile_inc, m

gen inc_low=inrange(quantile_inc,1,1)
replace inc_low=. if missing(quantile_inc)

gen inc_high=inrange(quantile_inc,2,2)
replace inc_high=. if missing(quantile_inc)

gen inc_all=1

keep county_fips inc_*

tempfile qinc
save "`qinc'"

restore

merge m:1 county_fips using "`qinc'", nogen
drop _merge

gen split_early=1 if effective_date>=mdy(3,17,2020) & effective_date<=mdy(3,26,2020) & effective_date!=.
gen split_late=1 if effective_date>= mdy(3,27,2020) & effective_date<= mdy(4,7,2020)
drop if state_fips==72
replace split_early=0 if split_late==1
replace split_late=0 if split_early==1

tab effective_date if split_early!=1 & split_late!=1, m

gen split_trumpe =1 if split_trump==1 & split_early==1
gen split_trumpl =1 if split_trump==1 & split_late==1
gen split_deme =1 if split_dem==1 & split_early==1
gen split_deml =1 if split_dem==1 & split_late==1



*Subset to merged observations from above datasets with non-missing counties and dates

keep if !missing(county_fips)
keep if !missing(date)


*Save copy of dataset for future analyses
save "${outdata}/countydata_foreventtime", replace




preserve
keep county_fips date_of_first_case effective_date
duplicates drop
format date_of_first_case %td
scatter effective_date date_of_first_case , scheme(s1mono) yti("Date of stay-at-home order") xti("Date of first case")
graph export "${figures}/firstcase_v_order.pdf", replace
restore 



*********************************************
***Event times based on stay at home order***
*********************************************

cd "${figures}"


local split "inc_all inc_high inc_low split_trump split_dem split_early split_late split_trumpe split_trumpl split_deme split_deml"

foreach i of local split {

local outcomes "distance visit pe_r_food pe_r_other pe_r_all log_pe_food log_pe_other log_pe_all"

foreach o of local outcomes {

tempfile `o'_`i'_s

preserve


*Subset to relevant sample and define label for graph

keep if `i'==1

sort county_fips date
*drop counties with no stay at home order 
drop if ever_stayathome==0

**keep from march 1
drop if date<mdy(3,1,2020)

*define the window for the eventtime graphs
local window 11
display `window'
local total (2*`window'+1)

**event time 0 is the day before effective date of order
gen eventtime=date-effective_date +1

replace eventtime = -11 if eventtime<-11
replace eventtime = 11 if eventtime>11
tab eventtime

**event time 10 is the effective date of order
**shift event time variables because stata doesn't like negative categories
replace eventtime=eventtime+`window'
tab eventtime

**omitted cat is -5
fvset base 6 eventtime 

gen beta=.
gen se=.
gen highci=.
gen lowci=.


reghdfe `o' i.eventtime deaths_daily_incr log_deaths cases_daily_incr log_cases [aweight=pop_estimate_2018], absorb(date county_fips days_since_first_case) cluster(state_fips date)



**pull the regression results for the figure
matrix b=e(b)
matrix V=e(V)

collapse (mean) baseline = `o' if eventtime==10 & e(sample)==1 [aweight=pop_estimate_2018]
local baseline = baseline[1]

**create a dataset from the regression results
**the way I coded the figures makes it a bit of a pain to plot multiple event times on the same graph,
**though I agree that we should  
clear
set obs 23
gen n=_n
gen baseline=.
gen beta=.
gen se=.
gen highci=.
gen lowci=.
forval x=1/23 {
replace beta=b[1,`x'] if n==(`x')
replace se=sqrt(V[`x',`x']) if n==(`x') 
replace highci=beta+(1.96*se) if n==(`x')
replace lowci= beta-(1.96*se) if n==(`x')
}
sort n
rename n date
replace date = date-12
keep if inrange(date,-10,10)
replace lowci=. if beta==0
replace highci=. if beta==0
replace baseline = `baseline'

rename beta b_`o'_`i'_s
rename high h_`o'_`i'_s
rename low l_`o'_`i'_s
rename baseline base_`o'_`i'_s
drop se

save ``o'_`i'_s', replace
restore
}
}


local distance_lab "% change from pre-pandemic levels   "
local visit_lab "% change from pre-pandemic levels   "
local pe_r_food_lab "Revenue per establishment ($)"
local pe_r_other_lab "Revenue per establishment ($)"
local pe_r_all_lab "Revenue per establishment ($)"
local log_pe_food_lab "Log revenue per establishment ($)"
local log_pe_other_lab "Log revenue per establishment ($)"
local log_pe_all_lab "Log revenue per establishment ($)"


**************************************************************
***Event times based on stay at home order: same regression***
**************************************************************

cd "${figures}"


local split " inc_high  split_dem split_early "

foreach i of local split {

local outcomes "distance visit pe_r_food pe_r_other pe_r_all log_pe_food log_pe_other log_pe_all"

foreach o of local outcomes {

preserve


sort county_fips date
*drop counties with no stay at home order 
drop if ever_stayathome==0

**keep from march 1
drop if date<mdy(3,1,2020)

*define the window for the eventtime graphs
local window 11
display `window'
local total (2*`window'+1)

**event time 0 is the day before effective date of order
gen eventtime=date-effective_date +1

replace eventtime = -11 if eventtime<-11
replace eventtime = 11 if eventtime>11
tab eventtime

**event time 10 is the effective date of order
**shift event time variables because stata doesn't like negative categories
replace eventtime=eventtime+`window'
tab eventtime

**omitted cat is -5
fvset base 6 eventtime 

reghdfe `o' i.eventtime##i.`i' deaths_daily_incr log_deaths cases_daily_incr log_cases [aweight=pop_estimate_2018], absorb(date county_fips days_since_first_case) cluster(state_fips date)

**pull the regression results for the figure
matrix b=e(b)
matrix V=e(V)

collapse (mean) baseline = `o' if eventtime==10 & e(sample)==1 [aweight=pop_estimate_2018]
local baseline = baseline[1]

**create a dataset from the regression results
clear
set obs 23
gen n=_n
gen baseline=.
gen beta1=.
gen se1=.
gen highci1=.
gen lowci1=.


***pulls in the main effects of the event time dummies (split=0)
forval x=1/23 {
replace beta1=b[1,`x'] if n==(`x')
replace se1=sqrt(V[`x',`x']) if n==(`x') 
replace highci1=beta1+(1.96*se) if n==(`x')
replace lowci1= beta1-(1.96*se) if n==(`x')
}

***adds main effect coefficients and the interactions (split=1)
gen beta2=.
gen se2=.
gen highci2=.
gen lowci2=.

forval x=0/22 {
 lincom `x'.eventtime + `x'.eventtime#1.`i' 
 replace beta2= r(estimate) if n==(`x'+1)
 replace se2 =   r(se)   if n==(`x'+1)
 replace highci2 =  r(ub)  if n==(`x'+1)
 replace lowci2 =  r(lb)  if n==(`x'+1)
}

 
sort n
rename n date
replace date = date-12
keep if inrange(date,-10,10)
replace lowci1=0 if beta1==0
replace highci1=0 if beta1==0
replace lowci2=0 if beta2==0
replace highci2=0 if beta2==0
replace baseline = `baseline'

if "`i'"=="inc_high" {
graph twoway ///
 (scatter beta2 date , scale(*1.4) c(l)  msize(vsmall) color(gs7) lpattern(shortdash)  ) ///
 (rarea highci2  lowci2 date ,  color(gs7%15) ///
  xline(0,lcolor(gs8)) yline(0,lcolor(gs8))  ytitle("``o'_lab'") ysc(titleg(0.1cm))  xtitle("Days relative to event") ///
  scheme(s1color) ti("", span linegap(0.01) size(medsmall))) ///
  (scatter beta1 date , c(l)  msize(vsmall) color(gs11)) ///
 (rarea highci1 lowci1 date ,  color(gs11%15) ///
 legend(order(1  "High income" 3 "Low income") region(lcolor(none)) ))  
 graph export "${figures}/`o'_incsplit_1reg.pdf", replace
}


if "`i'"=="split_dem" {

graph twoway ///
 (scatter beta1 date , scale(*1.4) c(l)  msize(vsmall) color(gs7) lpattern(shortdash)  ) ///
 (rarea highci1 lowci1 date ,  color(gs7%15) xline(0,lcolor(gs8)) yline(0,lcolor(gs8)) ///
 ytitle("``o'_lab'"  ) ysc(titleg(0.2cm))  xtitle("Days relative to event") scheme(s1color) ti("", span linegap(0.01) size(medsmall))) ///
  (scatter beta2 date , c(l)  msize(vsmall) color(gs11)   ) ///
 (rarea highci2 lowci2 date ,  color(gs11%15) ///
 legend(order( 1 "Republican (2016)" 3 "Democrat (2016)") region(lcolor(none)) )) 
 graph export "${figures}/`o'_trumpsplit_1reg.pdf", replace
}

if "`i'"=="split_early" {

graph twoway ///
 (scatter beta2 date , scale(*1.4) c(l)  msize(vsmall) color(gs7) lpattern(shortdash)  ) ///
 (rarea highci2 lowci2  date ,  color(gs7%15) xline(0,lcolor(gs8)) yline(0,lcolor(gs8)) ///
 ytitle("``o'_lab'" ) ysc(titleg(0.1cm))  xtitle("Days relative to event") scheme(s1color) ti("", span linegap(0.01) size(medsmall))) ///
  (scatter beta1 date , c(l)  msize(vsmall) color(gs11)) ///
 (rarea highci1 lowci1 date ,  color(gs11%15) ///
 legend(order( 1 "Early" 3 "Late") region(lcolor(none)) )) 
 graph export "${figures}/`o'_timingsplit_1reg.pdf", replace
 
}

restore
}
}

********************************
***Event times based on first case/death***
*********************************************************************

*Focus on data after February 24
drop if date<=mdy(2,23,2020)
local outcomes "distance visit pe_r_food pe_r_other pe_r_all log_pe_food log_pe_other log_pe_all"

*drop counties with no stay at home order 
drop if ever_stayathome==0
local medvar "case death"

**creates a variable with the first case or death
foreach m of local medvar {
gen first`m'_t=.
sort county_fips
bysort county_fips: replace first`m'_t=1 if `m's[_n]!=0 & `m's[_n-1]==0
replace first`m'_t=date*first`m'_t
bysort county_fips: egen first`m'_tt=min(first`m'_t)
bysort county_fips: egen first`m'=mean(first`m'_tt)
format first`m' %td
}
**drop counties where we don't have a 10 day period of no cases in our data 
drop if firstcase<mdy(3,7,2020)
local split "inc_all "

foreach i of local split {

foreach m of local medvar {
foreach o of local outcomes {
preserve
keep if `i'==1


tempfile `o'_`i'_`m'

*define the window for the eventtime graphs
local window 11
display `window'
local total (2*`window'+1)

**event time 0 is the day before first case/death

gen eventtime=date-first`m'+1

replace eventtime = -11 if eventtime<-11
replace eventtime = 11 if eventtime>11
tab eventtime, m

**event time 10 is the first case/death
replace eventtime=eventtime+`window'
tab eventtime

**omitted cat is -5
fvset base 6 eventtime 

gen beta=.
gen se=.
gen highci=.
gen lowci=.



reghdfe `o' i.eventtime deaths_daily_incr log_deaths cases_daily_incr log_cases [aweight=pop_estimate_2018], absorb(date county_fips days_since_first_case) cluster(state_fips date)


matrix b=e(b)
matrix V=e(V)

collapse (mean) baseline = `o' if eventtime==10 & e(sample)==1 [aweight=pop_estimate_2018]
local baseline = baseline[1]

**create a dataset from the regression results
clear
set obs 23
gen n=_n
gen baseline=.
gen beta=.
gen se=.
gen highci=.
gen lowci=.
forval x=1/23 {
replace beta=b[1,`x'] if n==(`x')
replace se=sqrt(V[`x',`x']) if n==(`x') 
replace highci=beta+(1.96*se) if n==(`x')
replace lowci= beta-(1.96*se) if n==(`x')
}
sort n
rename n date
replace date = date-12
keep if inrange(date,-10,10)
replace lowci=. if beta==0
replace highci=. if beta==0
replace baseline = `baseline'

rename beta b_`o'_`i'_`m'
rename high h_`o'_`i'_`m'
rename low l_`o'_`i'_`m'
rename baseline base_`o'_`i'_`m'
drop se
save ``o'_`i'_`m'', replace
restore
}
}
}




***Figure: for each outcome, plot the splits separately for stay at home event


local outcomes "distance visit pe_r_food pe_r_other pe_r_all log_pe_food log_pe_other log_pe_all"
foreach o of local outcomes { 

use ``o'_inc_all_s', clear

merge 1:1 date using ``o'_inc_high_s', nogen
merge 1:1 date using ``o'_inc_low_s', nogen
merge 1:1 date using ``o'_split_trump_s', nogen
merge 1:1 date using ``o'_split_dem_s', nogen
merge 1:1 date using ``o'_split_early_s', nogen
merge 1:1 date using ``o'_split_late_s', nogen

merge 1:1 date using ``o'_split_trumpe_s', nogen
merge 1:1 date using ``o'_split_trumpl_s', nogen
merge 1:1 date using ``o'_split_deme_s', nogen
merge 1:1 date using ``o'_split_deml_s', nogen

sort date

**makes confidence interval meet at zero for omitted category
**alternative: just have it not plotted at all at 0
replace h_`o'_inc_high_s =0 if date==-5
replace l_`o'_inc_high_s =0 if date==-5
replace h_`o'_inc_low_s =0 if date==-5
replace l_`o'_inc_low_s =0 if date==-5
replace h_`o'_split_trump_s =0 if date==-5
replace l_`o'_split_trump_s =0 if date==-5
replace h_`o'_split_dem_s =0 if date==-5
replace l_`o'_split_dem_s =0 if date==-5
replace h_`o'_split_early_s =0 if date==-5
replace l_`o'_split_early_s =0 if date==-5
replace h_`o'_split_late_s =0 if date==-5
replace l_`o'_split_late_s =0 if date==-5

replace h_`o'_split_trumpe_s =0 if date==-5
replace l_`o'_split_trumpe_s =0 if date==-5
replace h_`o'_split_deme_s =0 if date==-5
replace l_`o'_split_deme_s =0 if date==-5
replace h_`o'_split_trumpl_s =0 if date==-5
replace l_`o'_split_trumpl_s =0 if date==-5
replace h_`o'_split_deml_s =0 if date==-5
replace l_`o'_split_deml_s =0 if date==-5


graph twoway ///
 (scatter b_`o'_inc_high_s date , scale(*1.4) c(l)  msize(vsmall) color(gs7) lpattern(shortdash)  ) ///
 (rarea h_`o'_inc_high_s l_`o'_inc_high_s date ,  color(gs7%15) ///
  xline(0,lcolor(gs8)) yline(0,lcolor(gs8))  ytitle("``o'_lab'")  ysc(titleg(0.2cm)) xtitle("Days relative to event") ///
  scheme(s1color) ti("", span linegap(0.01) size(medsmall))) ///
  (scatter b_`o'_inc_low_s date , c(l)  msize(vsmall) color(gs11)) ///
 (rarea h_`o'_inc_low_s l_`o'_inc_low_s date ,  color(gs11%15) ///
 legend(order(1  "High income" 3 "Low income") region(lcolor(none)) ))  
 graph export "${figures}/`o'_incsplit.pdf", replace
 
graph twoway ///
 (scatter b_`o'_split_trump_s date , scale(*1.4) c(l)  msize(vsmall) color(gs7) lpattern(shortdash)  ) ///
 (rarea h_`o'_split_trump_s l_`o'_split_trump_s date ,  color(gs7%15) xline(0,lcolor(gs8)) yline(0,lcolor(gs8)) ///
 ytitle("``o'_lab'")  ysc(titleg(0.2cm)) xtitle("Days relative to event") scheme(s1color) ti("", span linegap(0.01) size(medsmall))) ///
  (scatter b_`o'_split_dem_s date , c(l)  msize(vsmall) color(gs11)) ///
 (rarea h_`o'_split_dem_s l_`o'_split_dem_s date ,  color(gs11%15) ///
 legend(order( 1 "Republican (2016)" 3 "Democrat (2016)") region(lcolor(none)) )) 
 graph export "${figures}/`o'_trumpsplit.pdf", replace
 
graph twoway ///
 (scatter b_`o'_split_early_s date , scale(*1.4) c(l)  msize(vsmall) color(gs7) lpattern(shortdash)  ) ///
 (rarea h_`o'_split_early_s l_`o'_split_early_s date ,  color(gs7%15) xline(0,lcolor(gs8)) yline(0,lcolor(gs8)) ///
 ytitle("``o'_lab'")  ysc(titleg(0.2cm)) xtitle("Days relative to event") scheme(s1color) ti("", span linegap(0.01) size(medsmall))) ///
  (scatter b_`o'_split_late_s date , c(l)  msize(vsmall) color(gs11)) ///
 (rarea h_`o'_split_late_s l_`o'_split_late_s date ,  color(gs11%15) ///
 legend(order( 1 "Early" 3 "Late") region(lcolor(none)) )) 
 graph export "${figures}/`o'_timingsplit.pdf", replace
 
 
 graph twoway ///
 (scatter b_`o'_split_trumpe_s date , scale(*1.4) c(l)  msize(vsmall) color(gs7) lpattern(shortdash)  ) ///
 (rarea h_`o'_split_trumpe_s l_`o'_split_trumpe_s date ,  color(gs7%15) xline(0,lcolor(gs8)) yline(0,lcolor(gs8)) ///
 ytitle("``o'_lab'")  ysc(titleg(0.2cm)) xtitle("Days relative to event") scheme(s1color) ti("", span linegap(0.01) size(medsmall))) ///
  (scatter b_`o'_split_trumpl_s date , c(l)  msize(vsmall) color(gs11)) ///
 (rarea h_`o'_split_trumpl_s l_`o'_split_trumpl_s date ,  color(gs11%15) ///
 legend(order( 1 "Republican: Early" 3 "Republican: Late") region(lcolor(none)) )) 
 graph export "${figures}/`o'_trumptimingsplit.pdf", replace
 
 graph twoway ///
 (scatter b_`o'_split_deme_s date , scale(*1.4) c(l)  msize(vsmall) color(gs7) lpattern(shortdash)  ) ///
 (rarea h_`o'_split_deme_s l_`o'_split_deme_s date ,  color(gs7%15) xline(0,lcolor(gs8)) yline(0,lcolor(gs8)) ///
 ytitle("``o'_lab'")  ysc(titleg(0.2cm))  xtitle("Days relative to event") scheme(s1color) ti("", span linegap(0.01) size(medsmall))) ///
  (scatter b_`o'_split_deml_s date , c(l)  msize(vsmall) color(gs11)) ///
 (rarea h_`o'_split_deml_s l_`o'_split_deml_s date ,  color(gs11%15) ///
 legend(order( 1 "Democrat: Early" 3 "Democrat: Late") region(lcolor(none)) )) 
 graph export "${figures}/`o'_demtimingsplit.pdf", replace
 
 export delimited  using  "${tables}/`o'_estimates_splits.csv", replace

 }


***Figure: for each outcome, plot the  different events on one figure (ignore split)
***Have to rescale events so they match
local outcomes "distance visit pe_r_food pe_r_other pe_r_all log_pe_food log_pe_other log_pe_all"
foreach o of local outcomes {
use ``o'_inc_all_s', clear
merge 1:1 date using ``o'_inc_all_case', nogen
merge 1:1 date using ``o'_inc_all_death', nogen

**makes confidence interval meet at zero for omitted category
**alternative: just have it not plotted at all at 0
replace h_`o'_inc_all_s =0 if date==-5
replace l_`o'_inc_all_s =0 if date==-5
replace h_`o'_inc_all_case =0 if date==-5
replace l_`o'_inc_all_case =0 if date==-5
replace h_`o'_inc_all_death =0 if date==-5
replace l_`o'_inc_all_death =0 if date==-5

**by death 
graph twoway ///
 (scatter b_`o'_inc_all_s date , scale(*1.4) c(l)  msize(vsmall) color(gs7) lpattern(shortdash)  ) ///
 (rarea h_`o'_inc_all_s l_`o'_inc_all_s date ,  color(gs7%15) xline(0,lcolor(gs8)) yline(0,lcolor(gs8))  ytitle("``o'_lab'") ysc(titleg(0.2cm)) xtitle("Days relative to event") scheme(s1color) ti("", span linegap(0.01) size(medsmall))) ///
  (scatter b_`o'_inc_all_death date , c(l)  msize(vsmall) color(gs11)) ///
 (rarea  h_`o'_inc_all_death l_`o'_inc_all_death date , c(l) clpat(tight_dot) msize(vsmall) color(gs11%15) ///
   legend(order(1 "Event: stay-at-home" 3 "Event: 1st death") row(1) region(lcolor(none)) )) 
 graph export "${figures}/`o'_byevent.pdf", replace
 
export  delimited  using "${tables}/`o'_estimates_byevent.csv", replace

}


clear
log close


